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Abstract 

We present the results of a finite-size analysis of the four dimen- 
sional abelian surface gauge model. This model is defined assigning 
abelian variables to the plaquettes of an hypercubical lattice, and is 
dual to the four dimensional Ising model. This last model is known to 
present a second order phase transition with mean field critical expo- 
nents. We have performed Monte Carlo simulations on several lattice 
sizes and high statistics. The analysis of the partition function zeroes 
and the specific heat scaling behaviour allowed us to estimate the crit- 
ical coupling P c as well as the critical exponents v and a. Our results 
are consistent with the second order critical exponents v = 1/2 and 
a = 0. The /3 C value is in perfect agreement with duality predictions 
from the 4-d Ising model. Nevertheless, the energy histograms show a 
seemingly non-vanishing double peak structure. The interface tension 
analysis suggests that this may be a finite size effect. 



2 Current address. 



Introduction 

The abelian surface gauge model was proposed some time ago in order 
to have a natural lattice translation of the three-index fields that appear in 
theories based on an antisymmetrical potential [1] (j)(x)^. Such theories are 
relevant in different domains as Chern- Simons studies, supergravity or the 
10-dimensional E s gauge theories [2]. Although different lattice strategies 
have been proposed [3], the surface gauge model seems to be a simple and 
numerically attainable way to study these theories. 

In a recent paper [4] a numerical analysis of the surface gauge model for 
the abelian Z(2) gauge group was performed for different space-time dimen- 
sionalities. In summary, the model is exactly solvable in three dimensions 
showing no phase transitions, evidences for a second order phase transition 
were found for the four dimensional case, whereas it was realized the exis- 
tence of a first order phase transition for higher dimensionalities, d=5 and 
6. Moreover, duality relations with several spin models were also studied. 
In particular, the four dimensional case is dual to the four dimensional Ising 
spin model, which is known to have a second order phase transition with 
mean field critical exponents. 

In this paper we present a high statistics finite-size analysis of the phase 
transition for d=4 on larger lattices. Our main purpose is to confirm the 
previous work, performed on rather small lattices, and to estimate the critical 
exponents of the theory. 

The model 

We denote by n a lattice site and by a link starting from point n in 
the direction \i of an 4-d hypercubical lattice. Starting from this point there 
are 4 positive links. A simple plaquette defined by directions /x, v is denoted 
by n^. We assign a gauge variable (o"i = ±1) to each plaquette. Elementary 
three-dimensional cubes are composed by six plaquettes {i £ d cu b e , i.e. the 
plaquette belongs to the perimeter of the cube). Three of them are the 
plaquettes associated to the combinations of the three index /x, v, r] related 
to the links starting from a given point n in the positive directions. The 
remainder three plaquettes are those that "close" the cube. The partition 
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function is denned, then, as: 

Z=^ e-? E ; E= £ E(cube); (1) 

{<7{} cube=l 

where 

£(cu&e) = l- J] ^ = ±1, (2) 

with /3 > being the inverse temperature in natural units. From now on, 
< . > will denote thermal averages of an observable in the canonical distri- 
bution. 



Duality 

A duality transformation [5] relates the four dimensional surface gauge 
model to the four dimensional Ising model. This last model has a second 
order phase transition with mean field critical exponents [6]. From series 
analysis [7], the critical temperature of the four dimensional Ising model has 
been estimated to be (3* = 0.149 65 ± 0.000 05. This value is in agreement 
with the numerical simulations of [6]. 

The duality transformation between these models implies the following 
relation between the critical couplings 

1 e 2/3 = + 1 

where (3* is the critical coupling of the four dimensional Ising model and j3 c 
is the surface gauge model one. 

From this equation, one can estimate the numerical value of the critical 
temperature for the surface gauge model giving j3 c = 0.953 44 ± 0.000 16, 
value in accordance with the numerical simulations of [4]. 

The finite size analysis of [6] found, for the 4-d Ising model, a second order 
phase transition with critical exponents a = 0.04±0.06 and v = 0.510±0.016, 
i.e. compatible with zero and with 1/2 respectively. These values for the 
critical exponents are, then, compatible with those expected from mean field. 
In his turn, duality implies also mean field critical exponents for the surface 
gauge case. Finally, it should be noticed that the Ising model has a local 
order parameter, the magnetization, that, due to the gauge invariance, is not 



2 



present in the gauge version. This reduces the number of critical exponents 
to be measured. 

Monte Carlo calculation 

Our numerical simulations with the surface gauge model have been per- 
formed applying the Metropolis algorithm [8] as it was done in Ref [4]. The 
Monte Carlo simulation took place for the following lattices and MC /3's: 



(L 


= 3; 


P = 


0.850), 


(L 


= 5; 


P = 


0.920), 


(L 


= 6; 


P = 


0.920, 0.930), 


(L 


= 9; 


P = 


0.940, 0.943), 


L = 


= 12; 


P = 


0.9474). 



To reach equilibrium, the first 50 000 sweeps were discarded from an 
initially ordered state, and further (from 500 000 to 1 000 000) sweeps were 
generated for measurements. After each sweep, the energy was measured and 
recorded in a time series file. Fig. la shows the first 50 000 energies recorded 
for L = 9 at (3 = 0.943. We checked that our results didn't depend on the 
initial configuration, studying how they changed when we discarded some of 
the bins into which we splitted the data. Besides, for L = 5 we repeated the 
simulation with a hot start. 

From the time series of E, it is straightforward to compute in the FSS 
region various quantities at nearby values of the Monte Carlo (MC) /3's by 
standard Ferrenberg-Swendsen reweighting [9]. To estimate the statistical 
errors, the time-series data was split into 20 bins, which were jack-knived 
[14] to decrease the bias in the analysis of the reweighted data. 

From the previous work in [4], and some new test runs, the above MC /3's 
were known to be close to the (lattice dependent) effective critical couplings, 
/3 C (L), of the system. We understand by effective critical beta the one at 
which the specific heat of a finite system has a peak. 

In the two cases (L = 6, 9) in which our first MC run was not close enough 
to the effective critical /3 C (L) of the system, we used the Ferrenberg-Swendsen 
reweighting procedure to determine it more precisely, and then launched a 
second MC run at this newly found effective critical beta. 
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In order to combine the information contained in the two MC runs (at 
different /3's) that took place for L = 6 and the two for L = 9, we used 
a patching procedure introduced in [10]. This patching procedure allows to 
estimate a single density of states from several MC runs by means of reweight- 
ing and least-square fitting of the energy histograms of the runs. Once the 
procedure has been applied, one can investigate the critical properties of the 
model in the vicinity of the phase transition. 

Fig. lb shows the energy histograms of our Monte Carlo simulations 
(L = 9; (3 = 0.94) and (L = 9; (3 = 0.943). Fig. lc shows the proba- 
bility distribution Pl(E;/3) at /3 = 0.9432 that we obtain from the energy 
histograms in Fig. lb by means of the patching procedure and the reweighting 
techniques already mentioned. 

Partition function zeroes 

From equations (l)-(2), we see that the energy of our model changes in 
steps of 4 within the interval < E < 8L 4 . Defining the new variable 

u = exp(-4/3), (5) 

the partition function becomes, then, a polynomial of degree 2L 4 + 1 in u. 
Using the Newton- Raphson method, as described in [11], we calculated the 
two closest zeroes to the real axis of the complex u-plane. The zeroes are 
given in Table 1. We have denoted the zeroes by u°, where the subindex i 
stands for the different zeroes, ordered according its distance to the real axis. 
For instance, u° corresponds to the zero with the smallest imaginary part. 

Let us define u c = u(/3 c ). It was first shown in [12] that, for L big enough, 
the finite-size scaling analysis of the partition function zeroes is given by 

u°(L) = u c + AL- 1 ^ [1 + 0(L-")}, cu>0, (6) 

where A is a complex constant. 

This result has been used both for first and second order phase transitions 
as well as spin systems and gauge theories [13, 10]. 

For the present case we use the imaginary part of the zeroes only. It turns 
out to be more effective. The reason is that the L dependence of the real 
part is rather weak and consequently has a worse signal-to-noise ratio than 
the imaginary part. This makes (6) less favorable, in this case, than just its 
imaginary part. 
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Table 1: Table of the 1 st and 2 nd zeroes of the partition function 



L 




Im{u\) 


Re(u° 2 ) 


Im(u° 2 ) 


3 


0.031630(34) 


0.008984(14) 


0.03143(12) 


0.020219(92) 


5 


0.025197(18) 


0.0022220(97) 


0.026316(72) 


0.004851(61) 


6 


0.024098(21) 


0.0014268(86) 


0.024905(43) 


0.003152(37) 


9 


0.022963(11) 


0.0005793(56) 


0.024205(68) 


0.001139(12) 


12 


0.022566(12) 


0.0003133(65) 


0.022789(42) 


0.000643(32) 



Due to the duality relations between our model and the 4d Ising model, 
we expect to obtain a result compatible with v = 0.5. 
In first approximation one can resort to 



Im{ul{L)) = Im(A)L- 1/v , 
which for pairs of lattices allows to estimate v{L, L') from 



v(L,L') 



In (L'/L) 



ln(Im(u%L'))/Im(u%L))) 



(7) 



(8) 



Our results are shown in Table 2. For increasing min(L, L') there seems 
to be a trend v{L, L') — > 0.5, as expected from duality. Our best estimate is 
(matching our largest lattices L = 9 and 12) 



v = 0.469(17), 



(9) 



which is consistent with 0.5 within two error bars. 

On the other hand, a linear regression fit of (7) with L = 6, 9, and 12 
gives us 

v = 0.4531(46) (10) 

with a goodness of the fit [11] of Q = 0.36. The lattices L = 3 and 5 worsen 
the fit. 
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Table 2: Estimates of v{L, L') for pairs of lattices of size L and L' 



L\L' 


5 


6 


9 


12 


3 


0.3657(12) 


0.3768(13) 


0.4008(14) 


0.4131(26) 


5 




0.4116(69) 


0.4372(34) 


0.4469(48) 


6 






0.4498(57) 


0.4572(65) 


9 








0.469(17) 



A first order phase transition would have a critical exponent v = 1/d = 
0.25 (d : spatial dimension). Our estimation of v, then, excludes the possi- 
bility of a first order transition. 

Specific heat 

As was mentioned before, we compute quantities at nearby values of the 
MC /3's by standard reweighting, and use jack-knived bins to take care of the 
statistical errors. In this fashion we determine the maxima of the specific 
heat 

C(/3) = ^(< E 2 > - < E > 2 ), (11) 

where N c is the number of three-dimensional cubes in lattice. 

The results are shown in Table 3, together with their effective critical /3's. 

Let us also define f3 zero (L) for the first zero of the partition function. 
From eq. (5) we have 

(3 zero (L) = - 1 - In [Re{ul{L)f + Im^L)) 2 ]. (12) 

It can be considered as another effective critical /3 C (L). Table 3 compares it, 
for different lattice sizes, with the specific heat effective critical /3's. 

For a second order phase transition, the locations of the specific heat (or 
any other set of effective critical /3's) are expected to scale as [15] 

= Pc + aL-V", (13) 
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Table 3: Table of the effective critical p zeT ° , the specific heat (C max ), and its 
-pseudo critical coupling 



L 


j^zero 


P c 

rmax 


n 

^max 


3 


0.85371(25) 


0.85277(25) 


1.7040(41) 


5 


0.91929(18) 


0.91825(17) 


2.514(18) 


6 


0.93097(22) 


0.92907(20) 


2.772(23) 


9 


0.94339(12) 


0.94300(11) 


3.106(40) 


12 


0.94780(13) 


0.94760(13) 


3.277(97) 



with a being a constant. If the temperature-driven transition of the model 
happened to be of first order, the scaling would have the same form, with v 
replaced by 1/d (d: spatial dimension). 

Assuming v = 1/2, we can use eq. (13) to obtain estimates of the critical 
coupling /3 c (oo) from linear fits in 1/L 2 . Our results from fits to the data 
(see Table 3) of the four and three largest lattices respectively are 

P c = 0.95390 ± 0.00012 (from (3° ax ) (14) 
P c = 0.95389 ± 0.00016 (from p zero ). (15) 

If instead of assuming v = 1/2 we use our best estimate v = 0.469 from eq. 
(9), the fits yield 

P c = 0.95294 ± 0.00012 (from p° ax ) (16) 
P c = 0.95268 ± 0.00016 (from p zero ). (17) 

The agreement between these results and the estimate based on the duality 
relation with the 4d Ising model, P = 0.95344(17), is slightly better when 
taking v = 1/2. 

We now turn to the specific heat which is usually the most difficult quan- 
tity to analyze. The reason is that the critical divergence is generally rather 
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weak and regular background terms become important. Recalling our v % 0.5 
and assuming hyperscaling to be valid for d = 4 as well, we expect 

a = 2-dv^0. (18) 

The corresponding FSS prediction is then [15] 

C max (L) = C((3Z ax (L),L) = B + B.lnL. (19) 

A linear fit with the three largest lattices gives B = 1.36(18) and B\ = 
0.790(94) with Q = 0.6. See Fig. 2. 

In should be remarked, however, that the confirmation of a = 0(log) is 
not really conclusive. Due to the small range over which C max varies we 
can fit the data also with a simple power-law Ansatz C max oc L a / V , yielding 
a/v = 0.265(31) with Q = 0.46 [from eq. (18) we have a/v = 2/v - d, 
and using our best estimate v = 0.469(17), results in a/v = 0.26(16)]. We 
also tried a non-linear three-parameter fit to the more reasonable Ansatz 
Cmax = K + biL a / v , but our data was not good, and the error bars turn out 
to be too large to draw any conclusion. 

Notice that in a first order phase transition we would expect a/v = d = 4. 

Interfacial free energy 

The Ferrenberg-Swendsen reweighting procedure [9] allows to estimate the 
canonical probability distribution Pl(E; /3) for a given lattice size and j3 from 
the energy histogram hist(E; /3mc) obtained in the Monte Carlo simulation 
(if our j3 of interest is close enough to the simulation /3). First order phase 
transitions are characterized by a double peak structure in the probability 
density Pl(E;/3) sufficiently close to the transition point. 

Fig. lc shows our estimated probability distribution for L = 9 at j3 = 

0. 9432. It has a typical first order transition double peak estructure. All 
our lattice sizes suffer the same effect. At the same time, we know from 
the previous sections that our estimates for the critical exponents v and a 
exclude a first order phase transition, and are consistent with a 4d Ising 
model like second order phase transition. 

In order to analyze this phenomena, we adopt the normalization P™ ax = 

1, and define an effective critical point (3^ by the requirement that both 
maxima are of equal height. If E L ' max < E L ' max , we have, then: 

L = f l{E l \P c ) = i J L = f l{E l ;P c ) = L. (ZU) 
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We introduce also P™ m = P L (E™ m ; (3%) as the minimum probability density 
between the two peaks. Fig. lc shows it for L = 9. 

Following [16], we argue that in order to be in the presence of a first 
order phase transition, P£ im has to decrease as L grows. More specifically, 
for a first order phase transition the interfacial free energy per unit area (or 
surface tension / s ) separating the coexisting ordered and disordered phases 
is given [17] by the limit L — > oo of the quantity 

2/£ = -^lln(Pr n )- (21) 

Table 4 presents our results for P£ im and /£ together with their corre- 
sponding flc- 

Table 4: Table of /3f , P£ m and the corresponding 2f s L 



L 




Dmin 
L 


2/£ 


3 


0.83870(22) 


0.4208(52) 


0.03206(46) 


5 


0.91885(22) 


0.609(32) 


0.00397(42) 


6 


0.92939(30) 


0.609(37) 


0.00230(28) 


9 


0.94320(10) 


0.557(24) 


0.000804(60) 


12 


0.94773(17) 


0.543(93) 


0.000355(69) 



Following [17, 18] we perform the FSS fit 

2f s = 2f' L + c/L*- 1 (22) 

with our four larger lattices in order to extract the L — > oo limit of eq. 21. 
It is depicted in Fig. 3 and gives 

2f s = 0.000103 ± 0.000070, (23) 

with c = —0.491 ± 0.045 for the constant, and a goodness of the fit Q = 0.8. 
Notice that this result is compatible with 2f s = within 2 error bars. In 
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fact, we have also tried the fit 

2fl = -c/L d -\ (24) 

and obtained c = -0.541 ± 0.029 with Q = 0.4. 

Based on these results, we argue that, in the infinite volume limit, the 
surface tension disappears, and we are left with a continuous phase transition 
in the critical point /3 C , as we expected from duality. 

Conclusions 

We have performed a fairly detailed Monte Carlo analysis of the 4-d 
abelian surface gauge model. Our estimates for the critical coupling (/3 C = 
0.95390(12), from the peak location of the specific heat, and (3 = 0.95389(16) 
from the closest zero to the real axis) agree with the estimate based on the 
duality relation with the 4d Ising model. 

From the finite-size scaling analysis of the first partition function zeroes 
we have extracted critical exponents v which discard a first order phase 
transition, and that manisfestly show a trend towards the v = 0.5 4d Ising 
result. Our best estimate is i/(9, 12) = 0.469(17). 

The analysis of the specific heat peaks is more problematic because our 
data is compatible with a critical exponent a = (as should be expected from 
duality), but also with a power-law scaling a ^ 0. In fact, we are unable to 
discriminate between the two Ansatz unless we go to larger lattices (L = 16 
would probably be enough). 

Altogether, considering our estimates for /3 C , v, and a, we can exclude a 
first order phase transition. It is arguable that they are a good Monte Carlo 
evidence in support of the duality relation between the surface gauge model 
and the 4d Ising model. 

The typical double peak structure that our probability densities present 
are somehow puzzling, because as we have seen a 4d Ising model like second 
order phase transition seems to hold for the model. A further finite-size 
analysis of the P™ m 's reveals that in the limit L — > oo the surface tension is 
compatible with zero, giving place to a continuous transition. 

We conclude, then, that the system experiments actually a second order 
phase transition with mean field critical exponents, as predicted by duality. 
The intriguing behaviour that originates the double peak structure seems 
to be due to an added phenomena that may be understood in terms of the 
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topological excitations of the model, i.e. the monopoles. In the QED case, 
these structures are 1-dimensional, and it has recently been noticed [19] that 
the imposition of periodic boundary conditions originate an unphysical sta- 
bility of these loops on the hypertorus, inducing a first order transition. In 
the present case, however, the topological excitations are point-like and this 
effect will not be present. Nevertheless, these excitations can manifest other 
phenomena, as condensation or site percolation, that can be responsible of 
the strange metastability over the phase transition. Note that this image is 
consistent with the fact than in five dimensions, where the topological ex- 
citations will become 1-dimensional, this model exhibits a clear first order 
phase transition. Work is in progress on this direction. 
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Figure Captions 



la Time history of the first 50 000 energy measurements at (3 = 0.943 on 
the L = 9 lattice. 

lb Energy histograms of our simulations (L = 9; j3 = 0.94) and (L = 
9; = 0.943). 

lc Canonical probability distribution Pl(E;(3), normalized by p^' max = 
p2,max ^ p _ Q_g^32_ Jt nas been estimated from the energy histograms 
in Fig. lb by means of the patching procedure introduced in [10] and 
the reweighting techniques of [9]. 

2 Finite-size scaling analysis of the specific-heat maxima (J max- Also 
shown are least-squares fits to a logarithmic Ansatz, C max = B + 
B x InL (with B = 1.36(18), B x = 0.790(94)), and to a pure power-law 
Ansatz, C max = cL a > v (with c = 1.73(10), a/u = 0.265(31)). 

3 Finite-size scaling fit 2, 2f s = 2f s L + c/L d ~ 1 (with c = -0.491(45)), 
and = 2fl + c\L d - x (with c = -0.541(29)), from the four larger 
lattices for the interfacial free energy. The dashed straight line is the 
asymptotic limit 2f s = 0.000103 obtained in the first case. 
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Fig. 1 can be requested via E-mail to baig@ifae.es or 
villanova@ifae.es 
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